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Introduction 


To understand the purpose of much of the mathematics in this unit, we need to 
look back to a group of units earlier in the course: Unit 6 (second-order 
differential equations), Unit 7 (oscillations) and Unit 8 (damping). In Unit 6 we 
studied the solution of differential equations, such as, for example, 


¥ + 4X + 10x = Scos 31, 


In Units 7 and 8 we studied topics in mechanics where such second-order 
differential equations are used. In Unit 7 we studied oscillations, and in particular 
simple harmonic motion, the ‘fabulous perfect’ oscillations described by equations 
of the form 


¥+ wx =0, (1) 
In Unit 8 we studied the model 
X + Qawk + wx =0 


which takes account of the damping exhibited in real mechanical systems. The 
behaviour of such a mechanical system under an externally applied sinusoidally 
varying force was also studied. This may be modelled by an equation such as 


X + 2awk + w?x = Acos(vt + ¢). (2) 


All the mechanical problems studied in Units 7 and 8 concerned the motion of a 
single object in one dimension. Suppose we wish to study the motion of a 
mechanical system involving two (or more) objects, or of a single object whose 
motion cannot be described by just one position coordinate. In such cases we may 
arrive (say by an application of Newton's laws) at a description of the motion of 
such an object involving linked differential equations in more than one variable, 
such as 

¥= -§) y 

bx + By ey 

j= —Gx - by. 
In the television programme associated with this unit we will examine the motion 
in a plane of an object whose position coordinates (x,y) are governed by these 
equations. We will also see how the equations can be solved. 


The substance of this unit is mathematics, although the motivation for much of its 
study comes from mechanics. The unit covers methods for solving certain systems 
of differential equations, such as (3) above, or 


2%, + 3%. = x, -4x. +e 

X, — 2%, =4x, + x, +sine, be 
In mechanics we are usually concerned with systems of equations, such as (3), 
involving second derivatives (and known as second-order systems). In this unit we 
will first study the simpler case of first-order systems (such as (4)) where only first 
derivatives occur. We will consider a method of solution for homogeneous first- 
order systems in Section 1, and a method for solving inhomogeneous first-order 
systems in Section 2. In Section 3 we study the solution of homogeneous second- 
order systems of a special form, similar to (3) above. Such systems may be 
regarded as an extension of simple harmonic motion (Equation (1)) to more than 
one variable. Finally, in Section 4, we look at systems of equations generalizing 
Equation (2) (which describes forced and damped motion) to more than one 
variable. We will not study the general solution of such systems, but will seek to 
extend to systems ideas such as ‘transient’ and ‘steady-state solution’ important in 
Unit 8. We shall also see how to compute steady-state solutions for such systems. 
The steady-state solution is often all that is required in mechanics. 


You have seen how systems of linear algebraic equations may conveniently be 
written using matrices. Systems of linear differential equations may also be 
expressed in matrix form. We shall find that the ideas of eigenvalue and 
eigenvector of a matrix, introduced in the preceding unit, are of value in solving 


We discuss the modelling of 
such mechanical systems in 
Unit 24. 


These linked differential 
equations were called 
‘coupled’ differential equations 
in Unit 15. 


Although the term ‘system’ is 
usually used in this context we 
could equally well talk of sets 
of differential equations. 


The term ‘homogeneous’ is 
defined in Subsection 1.1, and 
has much the same meaning 
as for a single differential 
equation. 
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systems of differential equations. Thus in this unit you will need to bring together 
ideas from earlier units on two apparently diverse topics: differential equations 
and matrices, 


Study guide 


The sections of this unit are not of equal length. You will probably find that 
Sections | and 2 take longer to study than Sections 3 and 4. 


There is no audio-tape for this unit. The television programme is part of Section 3, 
and the reading specifically associated with the programme is in Subsection 3.1. It 
is not necessary to have studied Sections 1 and 2 before watching the programme; 
however, you should read the pre-programme notes in Subsection 3.1 before 
viewing the programme. You will also find it helpful (but not essential) to have 
read the introduction to this unit and Subsection 1.1 before viewing. 


When studying this unit you will need to be familiar with the following material 
from previous units: 


Unit 2, Subsection 4.2 (the integrating factor technique); this is required for 
Subsection 2.1 of this unit. 


Unit 21, particularly Sections 1 and 2 (eigenvalues and eigenvectors); this is 
required throughout this unit. 


Unit 6, Section 4, or Unit 8; this is required for Section 4 of this unit, where we 
need: 


(i) the structure of the solution of a forced oscillation equation of the form (2) 
above; in particular, the ideas of ‘transient’ and ‘steady-state solution’; 


(ii) the calculation of steady-state solutions by the method of phasors. 


Unit 2, Section 2—reviewed in Unit 19, Section 1—(Euler’s method of numerical 
approximation); this is required for Subsection 2.2 of this unit. 


1 First-order systems 


1.1 Matrix notation and terminology 


As I mentioned in the introduction, a system of linear algebraic equations may 
also be written in matrix notation. For example 


3x;-— x2+2x3= 5 a | 2) |x } 
X, — 7x2 + 4x3 = -1 may be written 1-7) 4//x2] = |-1 
X2+ X3= 0 0 1 1} | x3 0 
or, more concisely, 
Ax=b 
3 =1 2 
where A is the matrix |1 —7 4), x is the vector [x, X2 x3]' and bis 
0 1 a 


the vector [5 —1 oy". 


Let us now consider the possibility of writing a system of differential equations in 
matrix form. 


The system 
24,4+3%)= x,-4nn +2 | 


2 . 7 (1) 
X, — 2X, =4x,+ x2. +sint 


6 MST204 22.1 


can also be written 


Lalli} 


To take full advantage of matrix notation we would like to write this as 


Ax = Bx + h(t), 
2 3 1 =4 2 
where A = ; ah B= ; fxs the column vector [x, x2)", and 


h(t)=[e _sint)’. Here, the entries in the column vector x are functions. We have 
extended the idea of differentiation to such vectors in a natural way by writing 


x= [i X2]". As far as vector algebra is concerned, column vectors of functions 
can be treated just like ordinary vectors. 
Example 1 


Suppose f(t) = (0? 3+4 eJPandg()=[2? 4  e% JF 
Then 
(f+ g) (0) = [3:7 T+t e@t+e*y? (addition) 


—-=(2 1 ef. (differentiation) 


Exercise er ) 
(i) Find Sam where g is defined in Example 1. 


(ii) If x = ae“, where a is a constant vector, find %. 
[Solution on p. 32} 


Using such vectors of functions, we can readily write suitable systems of 
differential equations in matrix form. 


Example 2 
Write 
¥, —3¥,=2x, -— x, 
28, + 4%. = x, + 3x, + sind } 
in matrix form. 
Solution 


In matrix form System (2) becomes 


AX = Bx + hit) 
2 =1 ; 
= j 3 x= fx, x2)" and 


Ree 1. -3 
where A = 
2 
The matrix notation need not be restricted to systems of two equations in two 


4 
h(t)= [0 — sin4r]’. 

unknowns, nor to matrices with constant entries. 
Example 3 


Write the following system of equations in matrix form. 
(X, — X3 = sine 


Xp +X) +x, =0 (3) 


(sin 1) 


Solution 


First write the system as 
1X, —x;=sint 


X2 + Xa + ty =0 


—X, + (sin rx, —2x; =0. 
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Then it is easily seen to be equivalent to 
A(x + B(t)x = h(t) 


t 0 0 i) 0 =1 
where A(t)= | 0 4 1),Bi)=|¢ O O|,x=[x, x2 x3)" 
—1 sint 0 0-2 «#0 


and h(t)= [sine 0 Ojf. 


However, not every system of differential equations can be written in matix form, 
in the way we have in these examples. For example the system 


XP 4X3 = xy + x2 
d } (4) 
Xk. =X, — x 
could not be. 
In this unit we shall be concerned almost exclusively with systems of the form 
Ay (0)X + Aa(t)x = h(r) (5) 
or of the form 
Ay (O¥ + A2(0)% + As(1)x = h(e) (6) 


(where A,(¢) is not the zero matrix in either case). Such systems are said to be 
linear. Equation (5) is a linear first-order system (it contains no derivative higher 
than first-order). Equation (6) is a linear second-order system. These definitions of 
‘linear’ and ‘order’ are extensions of those for a single differential equation. We 
shall usually be concerned with constant-coefficient systems, in which all the entries 
in the matrices A,(t), A2(t) and A3(t) are constants (though h(t) may depend on 1). 
If h(t) = 0 in (5) or (6), then we have a homogeneous linear system. 

Exercise 2 

Which of the systems (1)-(4) in the text above are linear? Of those that are linear, say 


which (if any) are homogeneous, and which (if any) are constant-coefficient, and give the 
order in each case. 


[Solution on p. 32] 
Exercise 3 
(i) Write the system 
X, + 2x, = 3%, — x2 
Xs =X) $x2 + x5 
%2 =k, +X5 
in matrix form, 
{ii)_Is this system constant-coefficient? Is it homogeneous? What is its order? 
[Solution on p. 32] 
In the remainder of this section we are concerned with methods of solution of 
linear, constant-coefficient, first-order systems of differential equations. In each 


case the first step in the method of solution is to write the system in normal form 
as defined below. 


A linear constant-coefficient first-order system is in normal form 
when it is written in the form 


& = Bx + k(t). 


The general constant-coefficient linear first-order system is 

Aix + Aox = h(i), 
This can be written in normal form so long as the matrix A, is non-singular. For 
then we have 


& = —Ay'Aax + Ay h(t), 


which is in normal form. 


‘Sce the introduction to Unit 6. 
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Exercise 4 
Where possible, write each of the systems of differential equations below in normal form. If 
it is not possible, explain why. 
(i) X + 2%. + 5x, =8 

2%, — X) + 10x, =2 
(ii) Xy +X. + Sx, =0 

2x, + 2%. +%,-22 =0 X, —%2 =x, + 2x, 
[Solution on p. 32] 


1.2 The solution of constant-coefficient homogeneous systems 


We now consider a method of solution for constant-coefficient, first-order, 
homogeneous systems. 


Method of solution 
Let us consider a constant-coefficient, first-order, homogeneous system in normal 
form: 

x = Bx. (7) 


We can derive quite a simple formula for the solution of (7), based on eigenvalues 
and eigenvectors. Let us first see where the formula comes from. To do this, we 
look for solutions of (7) which have the form 


x= ae“ (8) 
where a is a non-zero constant vector (and 4 is a constant scalar). fa =0we get x = 0, which is 
A a solution of (7), but not a 
If x = ae“, as in (8), then very interesting one. 


dae* (see Exercise 1(ii)). 
Thus (8) is a solution of (7), so long as 
jae" = Bae*. 
That is, so long as 
ja = Ba (where a # 0). (9) 


But this equation is just the definition of the eigenvalues and eigenvectors of the 
matrix B, So we see that x = ae“ is a non-zero solution of X = Bx so long as 4 is 
an eigenvalue of the matrix B and a is a corresponding eigenvector. 


Now the matrix B will (usually) have more than one eigenvalue, so x = Bx will 
(usually) have several solutions of the form x = ae“. When we have a 
homogeneous linear system of differential equations it is quite easily proved that 
any linear combination of solutions is again a solution. So if 4, and 4; are 
eigenvalues of B, and a, and a; are corresponding eigenvectors, then 

x = C,aye*" + Cra,e7" 
is also a solution of (7) (where C, and C, are arbitrary constants). 


We have found a number of solutions of (7). The only question that remains is 
whether we can construct the general solution of (7) by taking linear combinations 
of exponentials in this way. Unfortunately the answer is ‘not always’. The matrix B 
must satisfy the conditions given in Theorem | below. Fortunately, many of the 
matrices which occur in practice do satisfy the conditions of the theorem. 


Theorem 1 


Let B be a square, n x n, matrix. Suppose that B has n linearly 
independent eigenvectors, say a,,€2,...,8,, Corresponding to the 
(not necessarily distinct) eigenvalues A,,A2,...,A,. Then the 
general solution of the system of differential equations 


Bx 


Cyaje*" + Craze +... + CG 


You may notice a similarity to 
2 the general solution of a 
where C,,C>,...,C, are arbitrary constants. homogeneous second-order 
differential equation, in Unit 6. 
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Theorem 1, then, provides a procedure for finding the general solution of a 
homogeneous linear first-order system, so long as the system satisfies the 
conditions of the theorem. 


Procedure 1.2: To solve a linear, constant-coefficient, homogeneous, 
first-order system 


For this procedure to work we need to be able to express the 
system in the normal form 


x = Bx 


where B is an n x n matrix with n linearly independent 
eigenvectors, 


1, Express the system of equations in the normal form x = Bx. 


2. Find the eigenvalues A,,2,...,4, of Band a corresponding 
set of linearly independent eigenvectors a,,a2,...,8,- 


3. Write down the general solution of % = Bx in the form 
x = Cya,e*" + Cza,e*" +... + Cae" 


where C;,C>,...,C, are arbitrary constants. 


Example 4 
Find the general solution of 
X, =X, + Xp 
Xz = 3x, -— Ba 
Solution 


The system is already in normal form x = Bx where 


[ ql 

B= . 

3 -1 

We must find the eigenvalues and eigenvectors of B. 


The eigenvalues are found by solving the characteristic equation det(B — 21) = 0. 


Now 
det(B a=["5? , 
Vl lie ames 


=(1—4)(-1-—4)-3 
=B-4, 
So the eigenvalues are 2 and —2. 


To find corresponding eigenvectors we solve 


b=—2 1 u 
[ 3 aalle]-° 
putting 4 equal to each eigenvalue in turn. 
Case 4 = 2: we obtain the equations 
—u+ v=0 
3u—3v=0 


which have solutions of the form u = k, v = k. So an eigenvector corresponding to 
the eigenvalue 2 is [1 1 


Case 4 = —2: we obtain the equations 
3u+v=0 (twice) 
so an eigenvector corresponding to the eigenvalue —2 is [1 -3y. 


Step 1 


Step 2 


The general eigenvector is 
{k —_ k]’, but we only need 
one eigenvector here. 
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Since the conditions of Procedure 1.2 are satisfied (B is 2 x 2 with 2 linearly 
independent eigenvectors) we can write down the general solution: 


x= cof te + of _s]e*. 


That is © 


=Cet -2 
x1 = Cye* + Coe The choice of eigenvectors is 


ee 43 not unique. Other choices lead 
¥a = Cye* — 3C,e°™ to different looking, but 
where C, and C, are arbitrary constants. equivalent, forms of the 


general solution. 


It is sensible to check solutions to differential equations. This is readily done here. 
With x, and x; as above, we have 


%, = 2C,e" — 2C,e-, 
while 

X + x2 = 2C,e* — 2C,e7*. 
So X, = x; + Xz, as required. 
Similarly 

X2 = 2Cye* + 6C,e-** 

= 3x, — x2. 

Exercise 5 
Find the general solution of 

%, = 2x, + 3x, 

2 = 2x, + xy 
[Solution on p. 32} 


De cee ns Nie eee ne ope reee See Unit 21, Subsection 1.4, 
For Procedure 1.2 to work the matrix B must have n linearly independent for a discussion of ‘ 
eigenvectors, It does not matter if there are less than n distinct eigenvalues pte pod Sones ponding to 
provided there are enough linearly independent eigenvectors corresponding to each sl i 
eigenvalue to give a total of n linearly independent eigenvectors. Let us see an 

example of this. 


Example 5 

Find the general solution of the system 
X, =" Sx, + 3x5 
Xz = 3x, + 2x2 + 3x3 


X3 = —6x, — 4x3, 


Solution 
The equations are in normal form, as required. We have x = Bx with Step 1 
5) 0 3 
B=| 3 2 3 
—6 oe =4 
We need to find the eigenvalues and eigenvectors of B. Step 2 
The eigenvalues are found by solving the characteristic equation det(B — AI) = 0. 
We have 
S—A 0 Ss) 0 S-i a 
det(B—AI)=| 3 2-4 3 =-|2-4 a 3 
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2—A 3 3 
=|1e Sa 3 - e-af : | 
0 -6 -4-i4 
= (2-4) ((5 — a)(—4 — 4) + 18) 
=(2-A)(#? -1-2) 
= -(4-2P (+0). 
So the eigenvalues are 2 (a repeated eigenvalue) and —1, 
To find corresponding eigenvectors we solve 
54 0 3 u 
3 2-4 3 ||ol=0 
-6 0 -4-4||w 
putting A equal to each eigenvalue in turn. 
Case 1 = —1 gives 
6u + 3w=0 
3u + 3v + 3w=0 
—6u — 3w = 0. 
Putting u = k, the first and third equations give w = —2k. The second equation 
a a v a an eigenvector corresponding to the eigenvalue —1 is 
Case i = 2 gives 
3u+3w=0 
3u+3w=0 
—6u — 6w =0. 


So, putting w = k we have u = —k. Now v can take any value / say. Thus, any 
vector of the form [—k | k)" is an eigenvector. In particular 

(0 1 OjFand[-1 0 1)" are linearly independent eigenvectors 
corresponding to the eigenvalue 2. 

In this case, even though there are only two distinct eigenvalues, there are 


sufficient linearly independent eigenvectors (three) for the conditions of Procedure 
1.2 to be satisfied. Thus the general solution of the given system is 


1 0 ced | 
x=Cy] Lle'+Cz}1]e*+C3] Oe. 
ee 0 1 


In terms of x,, x2 and x;, the solution is 
x, = Cye™' — C3e”" 
X2 = Cie! + Cze* 
X3 = —2Cye™' + Cye™. 


Although the above example requires considerable calculation, there is little that is 
new in it. Most of it is a use of matrix methods you have met in earlier units. 


An example with complex eigenvectors 
In attempting to apply Theorem 1, it is quite possible that we may find that some 
or all of the eigenvalues and eigenvectors of the matrix B are complex. 
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Example 6 
Find the general solution of 


(10) 


Solution 


We have x = Bx where 


o ot 
B= . 
[4] 
To find the eigenvalues of B, look at the characteristic equation: 
= | 
-1 -d 
ie. 474+ 1=0. 


So the eigenvalues are i and —i. 


To find corresponding eigenvectors we solve 


[os all> 


putting 2 equal to each eigenvalue in turn. 
Case 2 = i gives 

-iu+v=0 

—u—iv=0. 


The first equation is i times the second and so both give v = iu. Thus an eigenvector 
corresponding to the eigenvalue i is [1 Cit 


Case 4 = —i gives 


iu+v=0 
—u+iv=0, 
The second equation is i times the first and so both give u = iv. Thus an eigenvector 
corresponding to the eigenvalue —i is [1 -i}’. 
So by Theorem 1, the general solution is given by linear combinations of x, and x; 
where 
1 ee cost + isint cost sint 
x=] Jet = = 5 = . [+i 
i ie —sint + icost —sint cost 
and 
ley ay" cost —isint cost] [sine 
X= _jet= Fer é é = 2 -ij l. 
—t —ie™ —sint — icost —sint COS f 
Since both x, and x, are linear combinations of [cost —sint)” = 4(x, +x) 


and [sint cost]? = paca — Xz), the general solution can equally well be written 


cost sint 
x-c[ & Jef i 
—sint cos ft 


For our purposes this is a much better way of writing the general solution because 
we are only interested in real solutions. 


The approach in the above example works because x, and x2 form a complex 
conjugate pair with real part [cost —sint]" and imaginary part 
(sint cost]. It turns out that systems for which B has complex eigenvalues 


Step 1 


Step 2 
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can always be treated in this way provided B is real (i.e. all the elements of B are 
real). The justification for this is provided by the following theorem. 


Theorem 2 


If the conditions of Theorem 1 are satisfied, and in addition B is 
real, then the terms in the general solution 


x = Cya,e*" + C.a,e" +... + C,a,e*" 
a is the vector obtained by 


are either real or occur in conjugate pairs ae“ and ae. taking the complex conjugate 
of each element of a. 


In Example 6 we'used Procedure 1.2 to find the general solution in the form 
x = C,ae" + Cyne" 


where a = [1 i)” and 4 = i. We then used the fact that the terms in this 
solution are conjugate to write the solution in the real form 


x = C,Re(ae“) + CzIm(ae*) 


where Re(ae“) = [cost —sinr]’ is the real part of ae“ and 
Im(ae“) = [sint —_cost]’ is the imaginary part of ae“. 


More generally, Theorem 2 states that if the general solution in Procedure 1.2 
contains complex terms then they occur in conjugate pairs. It is therefore possible 
to obtain a real form for the general solution by using the following procedure. 


Procedure 1.2(a): Applicable when the matrix B in Procedure 1.2 
real but has some complex eigenvalues 
Suppose that a is a complex eigenvector corresponding to the 


eigenvalue 7. Then replace the terms ae“ and He“ appearing in 
Procedure 1.2 by Re(ae“) and Im(ae“). 


The method of solution of x = Bx given in this section only works if B has enough 

linearly independent eigenvectors. In applications the matrix B often turns out to The eigenvectors of symmetric 
be symmetric. When this happens Procedure 1.2 is always applicable, for matrices were ciscussed relly 
symmetric matrices always have a full set of linearly independent eigenvectors. st the/ead of Unit 21, Section 1. 


Exercise 6 
Find the general solution of each of the systems below. 


(i) = 5x, + 4x, 
%= x 


(ii) 2%, = Sx, — 6x3 — 6x5 
X= —x, + 4x2 + 2x3 


3xy — 6x3 — dry 
#52 + 84-4 = (4 -1)(4 - 2)*). 
(iii) X, = —2x, + 2x, 
w= —M 


[Solution on p. 32] 


Summary of Section 1 


We can use matrix notation to represent suitable systems of differential equations. 
A linear first-order system is in normal form if it is written 


X=Bx+k. 


We can use Procedure 1.2 to solve a homogeneous, constant-coefficient system in 
the normal form, x = Bx, so long as B has enough linearly independent 
eigenvectors. We saw also in Procedure 1.2(a) how to write the solution in an 
explicitly real form, when some of the eigenvalues and eigenvectors of B are 
complex, 
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2 Further methods for linear first-order 
systems 


2.1 Constant-coefficient inhomogeneous systems 

Suppose we wish to solve an inhomogeneous system x = Bx + h(t), such as 
%,=2x,+3x,+ ¢ 

X2= 2x, + pare 0) 


So long as B has ‘enough’ eigenvectors (as in Theorem 1 of Section 1), it is 
possible to transform such a system into a much simpler system. We shall see 
shortly that by setting x; = 3y, + y2 and x; = 2y, — y2, System (1) becomes 


di=4nt+e \ 
Ja = —y2 — 20. 
Each equation in (2) can be solved by the integrating factor method of Unit 2, 


since each equation now involves only one unknown variable. We can then obtain 
the solution of (1) from that for (2). 


(2) 


Before describing the general procedure, let us see an example of how it works. I 
shall make a note of the steps involved in the margin. 


Example 1 
Solve System (1) above. 
Method 
The system may be written 
X = Bx + h(t) (3) 


3 


2 
where B = [; t 


|ana h(t)=[t 9 4r). 


The first step is to find the eigenvalues and eigenvectors of the matrix B. In fact we 
have already done this in the solution to Exercise 5 of Section | where we found 
that: 


(3 2)" is an eigenvector corresponding to the cigenvalue 4; 
ul —1)}" is an eigenvector corresponding to the eigenvalue —1. 


The next step is to form a matrix P, whose columns are the eigenvectors. So here, 


-[; i} 


We then know (by Theorem 5 of Unit 2/, Section 3) that P is non-singular and 
that 


rw 


That is, P~'BP is a diagonal matrix with diagonal entries equal to the eigenvalues, 
We now substitute x = Py in (3), where y is a new vector of variables. This gives 


Py = BPy + h(t). 


(Here we have used the fact that since all the entries in P are constant, 
differentiating Py with respect to ¢ gives Py.) 


Since P is non-singular we can multiply by P~! to obtain 
y=P'BPy + P"'h(e). 


Because P is a 2 x 2 matrix we can easily calculate its inverse: 


=1__1f-1 =1 
ia | ey VA 


Are there enou; 
eigenvectors? If so, form P so 
that P~'BP is diagonal. 


Put x = Py 


Calculate P~'. 
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Thus, in terms of the y variables, System (1) becomes 


Pll bei alle 


That is, 
Ue od 
Jo = —Y2 — 20. 


We now have a system that we can solve, because each equation can be solved by 
the integrating factor method of Unit 2. For example, the first equation can be 
written 

Fy —4y = 6 


Multiplying both sides by the integrating factor e~/*“ = e~* gives 


~4 


eo“ — 4y,) = te 
that is 
41) = te, 


Bye =te 
ae = 


Integrating both sides we obtain 
ye w= fear +, 


=(-}-teen 4 + Cy (integrating by parts), 
thus 
yi = Cre — dt — te. 
The second equation may be solved similarly to give (omitting the details) 
Yo = Cze™' — 28+ 2, 
Now the solution of the original system is x = Py. That is, 
[|-[3 | aie ai 
x2 2 -1)L Cre" -2r+2] 
Thus 
x, = 3Cye“ + Cre" — Yr + 
X2 = 2Cye" — Cre +94 -— 
So we have solved the original system (1). 


(The result of such a calculation should of course be checked, though I will not 
give details here.) 


We could have used a method similar to the above for the homogeneous systems 
in Subsection 1.2, but there the above method would have been unnecessarily 
complex. Here we need to make the transformation x = Py to get the system into 
a ‘diagonal’ form, that can be solved equation-by-equation using the integrating 
factor method. 


Solve the resulting equations 
mn 


y=P'BPy+P-th 
individually, by the integrating 
factor method, 


Now find x, from x = Py. 


Check 
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The procedure for inhomogeneous systems is set out below. 


Procedure 2.1: To solve a linear, constant-coefficient, 
inhomogeneous, first-order system 


For this procedure to work we need to be able to express the 
system in the normal form 


X = Bx + h(t) (3) 


where B is an n x n matrix with n linearly independent 
eigenvectors. 


1. (i) 
(ii) Form a matrix P whose columns are these eigenvectors; 
(iii) Find P~!, 

2. Put 


Find n linearly independent eigenvectors of B; 


x= Py 
in (3), to obtain 

y=Po'BPy + P~*h(1), (4) 
where P~'BP is diagonal. 


3. Find y by solving the equations in (4) individually, using the 
integrating factor method. 


4. Transform the solution for y back to a solution for x, using 


x= Py, 


Exercise 1 
Solve the system 
X= 2x, + 3x, + e* 
X_ = 2x, + x2 + 4e*, 
[Solution on p. 34] 


Exercise 2 
Solve the system 
%, = —2x, + 10x) + 48 
%, = 10x, — 2x, + 91e. 
[Solution on p. 34) 


2.2. Numerical methods 


We saw in Units 2 and 6 that differential equations cannot always be solved by 
analytical methods. In such a situation we may resort to obtaining numerical 
approximations to the solution. Techniques for doing this were discussed 
extensively in Unit 19. 


We are now faced with a similar situation for systems of differential equations. The 
analytical techniques described in this unit are only applicable to certain special 
types of systems and so the need to have numerical techniques for systems of 
differential equations is as pressing as it was for single differential equations. 


Although, for reasons of space, I cannot give a comprehensive treatment of the 
topic, it is useful to indicate here how we can go about constructing numerical 
schemes to generate approximate solutions to first-order systems of differential 
equations. We will look at a numerical method based on Euler’s method for a 


Euler's method was discussed 
in Units 2 and 19. 
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single differential equation. Recall that in Euler's method for a single differential 
equation we approximate the solution of the initial condition problem 


dx 
wae in(x, t) X = Xo when t = fo 


by starting with Xo = xp and using the recurrence relation 
Xp41 =X, + hm(X,,6,) 

to generate approximations X, to the true solution x(t,) at f, = to + rh. 

We can easily extend Euler’s method to a system of two equations of the form 
¥ 1, = my (X1,X2,0) 
Xz = Mz (X1,X2,t). 


We do this by using Euler's method to generate approximate solutions to each of 
these equations. Looking at the first, Euler's method gives an approximation X,, 
to x, (t,) where 


Kiyo = Xiy + hm (X17, X29 tr). 


Here X;, is the approximation to x,(t,) where (applying Euler's method to the 
second equation) 


Xaye1 = X2, + hm2(X,,,X2, 


This may seem a little complicated but in practice it is quite straightforward. Let 
us look at an example. 


Example 2 
Write down the recurrence relations corresponding to an application of Euler's 
method to the system 


¥, = tx, — Px, 
%_ = 3x, — cost, 
where x, = 1, x2 = —3atr=0. 
Solution 
Applying Euler's method to each equation in turn gives 
Xaeer = Xiy + WGK, — X27) 
Xap41 = X2, + h(3X,, —cost,). 
We have also that X\9 = 1, X2,9 = —3, and that ¢, = rh. 


(5) 


The recurrence relations (5) obtained in the Example can be used to generate 
approximations to the solution of the given system, step-by-step. Starting from the 
known values of X, 9 and X>\9, we use (5) to calculate X,,, and X2,. We then 
use (5) again to calculate X,,) and X2.2; and so on. The results obtained depend 
on the step-length h chosen, and will usually be more accurate the smaller h is. 


Example 3 
Use the recurrence relations (5) to calculate X,,. and X22 with h = 0.1. What 
values of the solutions do these approximate? 


Solution 
Putting r = 0 in (5), we obtain 
X11 =X1,0 + Ol (foX1,0 — 3X 2,0) 
X21 = X29 + 0.1 (3X19 — costo). 
Substituting X, 9 = 1, X29 = —3, and to = 0, we obtain 
Xi,=1 
X21. = —3 + 0.1(3 — 1) = -28. 
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Putting r = 1 1 (5) gives 
Xi2=X11 +01 0.1X,, —0.01X3,,) 
=1+0.1 (0.1 + 0.028) 
= 1.0128. 
X22 = —28 + 0.1 (3 —cos0.1) 
—2,5995. 


Since we are working with step length h = 0.1, X,,2 is the calculated 
approximation to x, (0.2), and X25 the approximation to x3(0.2). 


i 


Calculations such as this are time-consuming and would normally be performed 
on a computer, so I shall do no more by hand here. 


It is a straightforward matter to extend the above method to systems containing 
more than two equations. We can consider systems of equations which can be 
written in the form 


Xp = My (X4,X3y.0-5 Xml) 


Xz = M2 (X1,X24.--4 Xml) 


Ky = My(X 14X25 -0 4 Xm Oe 


If we introduce the vector notation x = [x; x ae x,]’ these equations 
can be written more concisely in the form 


X = m(x, f). 
The procedure for applying Euler's method to such systems is given below. 


Procedure 2.2: Euler’s method for systems of equations 
This procedure applies to systems of the form 

X= m(x, 1) 
with initial condition x = Xo at f = fo. 


1. Replace each equation in the system by a recurrence relation, 
using Euler's method. Thus the equation 


Xj = m,(x,t) 
is replaced by 
Xires = Xi + hm (X,,6,). (i=1,2,....) (8) 


Here X,=[X1, Xay «a:  Xa,fo inthe 
approximation to x(t,) where t, =f + rh. 


2. (i) Set Xo =xo. 
(ii) Choose a step length h. 


(iii) Use the recurrence relations (8) to calculate X,, then use 
(8) again to calculate X>, and so on. 


Exercise 3 


(i) Write down, in a form suitable for generating solutions, the recurrence relations 
corresponding to an application of Euler's method to 


X, = 3x, +4 
X2 =X — x2 -e, 


where x, = 5, x2 =2atr=0. 
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(ii) Use the recurrence relations you obtain, with step-length 0.1, to calculate 
approximations to x, (0.2) and x,(0.2). 


[Solution on p. 34) 


We have seen in previous units that care must be taken in the use of numerical 
methods. There is no time here for a discussion of the important topic of stability, 
which is essential to a proper understanding of numerical schemes such as those 
constructed above. Approaches other than those based on Euler’s method are 
possible; in fact we could develop methods for first-order systems based on any of 
the methods for a single first-order equation described in Unit 19, but, again, 1 will 
not pursue this topic here. 


Summary of Section 2 


To solve an inhomogeneous constant-coefficient system of differential equations 
x = Bx + h(t) where B is an n x n matrix with n linearly independent eigenvectors, 
we use Procedure 2.1. 


A method of constructing recurrence relations to generate approximate solutions 
to systems of differential equations was described. This is based on Euler's method, 
and consists of applying the approximation specified by Euler’s method to each 
equation in the system in turn. 


3 Second-order homogeneous systems of the 
form X = Bx: simple harmonic motion 


We now turn our attention to systems of differential equations which involve 
second derivatives. We shall be concerned solely with linear, constant-coefficient, 
second-order systems. Such systems arise frequently in modelling, especially in 
mechanics. A particularly important class of second-order systems are those that 
can be written in the form ¥ = Bx, and we study these in this section. 


Subsection 3.1 summarizes material covered in the television programme ‘Two-way 
motion’. The programme looks at the solution of a particular system of the form 

% = Bx and illustrates the sort of real situation in which such a system of 
differential equations may arise. In Subsection 3.2 we look at a procedure for 
solving systems of the form % = Bx. 


3.1 ‘Two-way motion’ (Television Subsection) 
Pre-programme notes 
In the programme we study the solution of the pair of differential equations 
—ix + by } 
3x — $y. 


This pair of differential equations provides an approximate description of the 
motion of a ball-bearing in a suitably shaped bowl (where x and y are the 
horizontal position coordinates of the ball-bearing). The programme consists of 
three parts: 


1. The actual motion of the ball-bearing in such a bowl; 


(ly 


2, A computer animation of the precise motion implied by Equations (1); 
3. The mathematical solution of Equations (1). 


(The computer animation enables us to see how accurately Equations (1) describe 
the actual motion of the ball-bearing. For various reasons the actual motion is 
only rather crudely approximated by Equations (1). For example the equations 
take no account of various factors tending to damp the motion.) 


For a discussion of stability 
see Unit 19. 


bs 


V 


Figure 1. We consider the 
motion of a ball-bearing in a 
bowl. The two horizontal 
components of acceleration, * 
and jj, are specified 
(approximately) by the 
Equations (1). 
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To solve the equations mathematically it is helpful to use vector notation. If we 
write r for the vector of variables [x y]' and A for the matrix 


[3 i] 
3-8) - 
then we see that the pair of differential equations (1) can be written 
¥=Ar. 
This equation is similar in form to the single differential equation 
5 = ky. 


For k < 0, this is the equation of simple harmonic motion. The similarity of the 
equations f = Ar and jj = ky is not simply a matter of mathematical formalism. 
For suitable matrices A, the vector equation ¥ = Ar does describe a generalization 
of simple harmonic motion to motion in more than one dimension. The modelling 
of mechanical systems displaying such motion is discussed in Unit 24 on normal 
modes. Our concern here is with the mathematical solution of the equation ¥ = Ar. 
However, in the programme we take an informal look at an example of two- 
dimensional motion—the ball in the bowl. This example is provided to help you 
to visualize the solution of the equation F = Ar which is obtained in the 
programme. 


The solution of the equation involves knowledge of the eigenvalues and 
eigenvectors of the matrix A. I will ask you to find these before viewing the 
programme. 


Exercise 1 
Calculate the eigenvalues and eigenvectors of 


a 
a-[4 ol 
[Solution on p. 35} 


Now view the television programme ‘Two-way motion’. 


Post-programme notes 
To solve a system of the form 
¥=Br (2) 
we may look for solutions of the form 
r=ae" (3) 
where a is a constant vector, and a ¥ 0. With r given by Equation (3), we have 
f=ypae and #=p'ae". 
Hence (3) is a solution of # = Br if, and only if, 
ae" = Bae“; 
ie. wa =Ba. 


With a # 0, this condition asserts that y? is an eigenvalue of the matrix B, and a 
is a corresponding eigenvector. 


If A is an eigenvalue of B and a is a corresponding eigenvector, we have two 
solution of ¥ = Br of the form r = ae”. These are 


r=ae“*'and r 


aens4" 


In the programme, we concentrated on the pair of Equations (1) on page 19, 
which are equivalent to f= Ar withr=[x —y]’ and 


“(td ° 


You will find this notation 
discussed further in Subsection 
11. 


Simple harmonic motion was 
discussed in Unit 7. 


The ball actually moves in 
three directions but its motion 
is completely determined by 
the behaviour of its two 
horizontal coordinates x and y. 


TV22 


We are not using exactly the 
same notation as in the 
programme. 
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The eigenvalues and eigenvectors of A are given in the solution of Exercise 1, 
where we found that: 


{1 —1]" is an eigenvector corresponding to the eigenvalue —1; 
tl 1]? is an eigenvector corresponding to the eigenvalue —$. 
So, the solutions of ¥ = Ar given by the eigenvalue —1 are 


More generally, a linear combination of these solutions is also a solution. That is 


te [i foe +d,e~"), 


where c, and d, are constants, is a solution. Use of Euler’s formula enables us to 
rewrite this as 


r= [_i}e cost + D, sint) (5) 


where C, and D, are arbitrary real constants. Similarly, corresponding to the 
eigenvalue —4, we have a solution 


v= [fem eet 
or 


r= [1 eos 4t + Dz sin 41). (6) 


Although (5) and (6) contain arbitrary constants, neither is the general solution of 
¥ = Ar. The general solution is a linear combination of both (5) and (6). That is 


r= mG cost + D, sint) + [1 }cssoet + D, sin 41). (7) 


The programme showed computer animations of the motion of a point in a plane, 
where the position coordinates (x, y) of the point satisfy the original pair of 
differential equations (1) (or, equivalently, # = Ar with A given by (4) and 

r=[x  y]’). 

Each of the Solutions (5) and (6) above corresponds to a particularly interesting 
motion of such a point. Looking at (5) we see that it gives the position 

r={x  y)" asa scalar multiple (C, cost + D, sint) of a fixed vector [1 —1]’. 
So, in the motion described by (5), the point moves along a straight line defined 
by the vector [1 —1]" (see Figure 2). The factor C, cost + D, sint varies 
sinusoidally, so that the point oscillates in simple harmonic motion along this line. 


Similarly Equation (6) implies simple harmonic motion along the line defined by 
the vector [1 1)". Notice also that Equations (5) and (6) show that the 
angular frequencies of these two special motions are in the ratio 1:3. 


The motions described by Equations (5) and (6) are only achieved under suitable 
initial conditions. The general motion of the point, under arbitrary initial 
conditions, is given by the general solution (7) of the equation F = Ar. But any 
solution of the form (7) is just the sum of a solution of the form (5) and one of the 
form (6). Thus any motion of the point can be seen as the sum of two special 
motions, one along each of the lines defined by [1 1J" and [1 —1)". The 
programme included a computer animation illustrating this fact. Figure 3 shows 
examples of the paths (orbits) that such a point may follow. The general motion 
described by the Equations (1), or their solution (7), is of the point tracing and 
retracing orbits of the form shown in Figure 3. 


Now consider the real bowl. By viewing the bowl from above we observe the 
behaviour of the horizontal x, y co-ordinates of the ball bearing. By comparing 


Euler's formula was discussed 
in Unit 5. It states that 


e” = cos0 + isin 0. 


y 


Figure 2. If the ball-bearing is 
released from either of the 

tions A or B, then the 

uations (1) imply simple 
harmonic oscillations along 
the straight lines /, or [, 
respectively. Such straight line 
motions were observed in the 
real bowl. 


(a) Motion of the point 
starting from rest at A. 


(b) Motion of the point 
starting from B in the 
direction shown. 


Figure 3. Examples of orbits 
determined by the Equations (1). 
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this view of the bow! with the computer animations we can judge how accurately 
Equations (1) describe the motion of the ball bearing. In the programme we found 
that the ball bearing could be made to describe oscillation along the lines defined 
by [1 1]’ and [1 —1]". However, the oscillations were damped rather than 
pure simple harmonic motion. Furthermore, the general orbits of the ball bearing 
in the real bowl were not of the form shown in Figure 3, for as well as damping, 
there were other inaccuracies in the real bowl. One was that the frequencies of the 
two special motions were not exactly in the ratio 1:3. We concluded the 
programme by trying to adjust the mathematical solution so as to make it give a 
more accurate description of the behaviour of the real ball-bearing. We did this by 
looking at computer animations obtained by adding simple harmonic motions in 
the same two directions [1 - 1]’ and{1 —1]" but with frequencies in a ratio 
which differed slightly from 1:3 (see Figure 4). 


Exercise 2 

An object moves so that its position, r = [x  y]", in a plane satisfies the equations 
¥= 4x + Sy 
P= hx — Wy, 


Give the directions (if any) in which such an object can describe simple harmonic motion in 
a straight line, and the frequency of any such motion. 


[Solution on p. 35) 


3.2 Solution of second-order systems of the form ¥ = Bx 
In the previous subsection we saw that 

x = ae" (4) 
is a solution of 

X = Bx (3) 


so long as yi? is an eigenvalue of B and a is a corresponding eigenvector. For each 
non-zero eigenvalue / of B, we have two solutions of (3): aev*’ and ae~v*", Aside 
from the fact that there are more terms in the solution, solving a second-order 
system of the form (3) is no harder than solving a first-order system. The general 
solution of (3) can usually be built up from solutions of the form (4). But as in 
Subsection 1.2 this method does not invariably work—it requires that the matrix 
B has enough eigenvectors. The following theorem and procedure for second-order 
systems of the form ¥ = Bx are very similar to those in Subsection 1.2 for first- 
order systems. 


Theorem 1 

Let B be an n x n matrix. Suppose that B has n linearly independent 
eigenvectors @;,2,...,@,, corresponding to the eigenvalues 

AyyAry ssn Fespectively. 


If all the eigenvalues are non-zero the general solution of the 
system of differential equations 


X = Bx 


x= ¥ a(cev*'+ de-v* (5) 


ret 
where c, and d, are arbitrary constants (1 <r <n). 


If any eigenvalue 4, is zero the general solution is obtained by 
replacing the rth term in (5) by (c, + d,t). 


In the example studied in the previous subsection the eigenvalues were negative 
and we had to express the general solution in terms of real functions. In general, 
the real form of the terms c,e**' + d,e~¥*"' in (5) depend on the nature of the 


Figure 4. By adjusting the 
ratio of the animation’s 
frequencies in the special 
directions, we can obtain a 
better correspondence between 
the animation and the motion 
in the real bowl. 
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eigenvalues. There are three cases to consider: (1) 4, real and positive, (2) A, real 
and negative and (3) 4, complex. (The case J, = 0 is considered separately in 
Theorem 1.) 
Case 1: A, real and positive 
In this case the square roots of 4, are real and the term in (5), 
a, (cer + de-V*9, 
is already explicitly real. 
Case 2: 4, real and negative 
Now the square roots of 2, can be written as 
if=%, and -i/=7,. 
Since both e’~** and e~'V-*" are combinations of cos J=het and sin J, I 
the term in (5) can in this case be written as 
a,(C,cos ,/—4, t + D,sin,/—3, 0). 


where C, and D, are combinations of ¢, and d, and are again arbitrary constants. 


Case 3: 4, complex 
In this case we can replace the complex exponentials in (5) in a way similar to that 
used in Procedure 1.2(a) for homogeneous first-order systems. 
We know from Theorem 1 that if A is a complex eigenvalue of B with 
corresponding eigenvector a then the functions 

ae¥*'and ae~¥* (6) 
will appear in the general solution (5). It turns out that if the matrix B is real then 
J will also be an eigenvalue of B, and 4 will be a corresponding eigenvector. This 
implies that the functions 

aev™' and ge7*7" (7) 
will also appear in the general solution. Now the functions (7) are complex 
conjugates of the functions (6) and so we can replace all four functions by the real 
functions 

Re(ae’*"), — Im(aev*"),—Re(ae~¥*),,_—and_—Im(ae~¥"), 
To summarize, we have the following procedure: 


Procedure 3.2: 
To solve a linear, second-order, system of the form ¥ = Bx 


For this procedure to work B must be an n x n matrix with n 
linearly independent eigenvectors. 


1, Find the eigenvalues 2,,2,...,4, of B and a corresponding 
set of linearly independent eigenvectors a,,a3,...,a,. 


2. Write down the general solution in real form. There are four 
cases to consider: 


(i) If all the eigenvalues are real and positive the general 
solution is 


LJ rr rc 
x= ¥ a(Ce'+ Dew, (8) 
m= 
(ii) If 4, = 0 replace the rth term in (8) by 
a,(C, + D,t). 
(iii) If A, is real and negative replace the rth term in (8) by 
a,(C,cos,/—4, t + D,sin/—2, 0). 


(iv) Ifan eigenvalue 2 is complex replace the functions ae*', 
ae~**" gee" and ge~ "in (8) by the functions 


Re(ae~*"), Im(ae**"), Re(ae~**') and Im(ae~¥*"), 
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The eigenvalues of the matrix B do not necessarily all fall into the same case, of 
course. The example below illustrates this. 


Example 1 

Find the general solution of the system 
Xy= Xy +3x2-— X3 
¥,=3x,+ 2+ X3 
X= 4x3. 

Solution 


We have a system of the form ¥ = Bx with 


1 3) 
B= |3 1 1}. 
oO) oO 4 


The first step is to find the eigenvalues and eigenvectors of B. 


The eigenvalues can be found by solving the characteristic equation 
det(B — AI) = 0. We have 


1-1 3 at 
det(@B-A)=| 3 1-A 1 


"| 
-Aa 
= (4—A)((l — AP - 3?) 
= -(44-A)(4—4)(2+2). 
So the eigenvalues are —2 and 4 (repeated). 


1-4 
=a-afs ; 


The eigenvectors corresponding to the eigenvalue —2 satisfy 
3 3 -I1\[u 0 
3 3 1} |v |=] 0}. 
oO 0 6) Lw. ) 


Thus w = 0 and u = —», and so an eigenvector is [1 —1 0)”. The 
eigenvectors corresponding to the eigenvalue 4 satisfy 


=—3| 3 =i} ) 
3 3 1} |e} =]0}. 
0 0 Oj}lw 0 


These equations reduce to 
—3u+3v-—w=0. 


Since we effectively have only one equation in three unknowns, we can find two 
linearly independent eigenvectors corresponding to the eigenvalue 4: for example 
fl 1 0]’ and [0 1 ai. 


By Procedure 3.2 the general solution is 
i 0 1 


x =| 1|(Cye + Die“) + | 1 |(Cze* + Dze~*) + |—1|(C3cos /2t + Dysin/2t) 
0. 3 0. 
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or alternatively 


x1 = Cye* + Dye“™ + Cycos/2t + Dy sin /2t 
X2 = (Ci + Ca)e* + (D, + Dz)e™™ — C3cos,/2t — Dysin\/2t 
x3 = 3(Cze" + Dye~'), 

Exercise 3 

Find the general solution of each of the systems below. 


2%, = — 12x; — 20x2. 
[Solution on p. 35] 


Summary of Section 3 


To solve a second-order, constant-coefficient, system of the form X = Bx, for 
suitable B we may proceed in a similar way to that described in Section 1 for a 
first-order system. The details are given in Procedure 3.2. 


Of particular significance in mechanics is the case where the eigenvectors of B are 
real and negative. In this case the model ¥ = Bx represents a generalization of 
simple harmonic motion (Unit 7), In the television programme we looked at an 
example of this: the two-dimensional motion of an object whose position vector 
r=[x  y]’ varied with time according to the system of differential equations 
¥ = Ar where the matrix A had negative eigenvalues. 


4 Forced oscillations 


4.1 The structure of the solution 


In the unit ‘Differential equations IP’ we introduced certain results about the 
solution of the forced oscillation equation 

¥ + 2xx + w?x = by cos vt + by sin vt. 
These results were used extensively in the mechanics unit on “Damped and forced 
vibrations’. 


In this section, we discuss the extension of these results to systems of differential 
equations. I will first remind you of the results for a single equation, by looking at 
an example. 


Example 1 
We consider the solution of the equation 
X¥ + 2% + Sx =2cos 31. (1) 


To solve this equation we first calculate the complementary function; that is, the 
general solution of the associated homogeneous equation 

¥ + 2x + Sx =0. 
I will not give details of this calculation. The complementary function is 

X_ = e-'(A.cos 2t + Bsin2t). 


Next we calculate a particular solution. We can do this by the method of phasors, 
in which we try x = Re(ze**) and calculate the complex number z. Since the 
phasor of 2 cos 3¢ is 2, we get here 


((3i)? + 231) + 5)z = 2, 


and so 
2 1 i 
i= a pe + 3i). 


Note the 2 on the left-hand 
side of part (ii), 


Unit 6, Section 4 


Unit 8 
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Thus a particular solution is 
Xp = —Fycos3t + yysin 30. 


The general solution of (1) is found by adding the particular solution and the 
complementary function, and is 


X=X.+Xp 
= e7'(A cos 2t + Bsin 2t) — fy cos 3t + Fysin 32. 
Now let us consider the behaviour of this solution for large values of t. The 


complementary function tends to zero for large t, because of the term e~'. Thus for 
large ¢ any solution of (1) is approximately equal to the particular solution 


Xp = —fycos 3t + yysin 3t. 


This is known as the steady-state solution of (1). We refer to the complementary 
function as a’ transient, since its effects die out as t becomes large. We saw in Unir 8 
that in many mechanical problems on forced vibrations, a knowledge of the 
steady-state solution is all that is needed, so long as we know that the 
complementary function is indeed a transient (ie. it tends to zero for large 1). 


In general, we have the following important information about the solution of the 
single equation 
4,X + aX + ax = by sinwt + bz cos wt. (2) 
(i) | The general solution of (2) is the sum of the complementary function and 
any one particular solution, 


(ii) So long:as each of a;, a; and as is positive, the complementary function is a 
transient (i.e. it tends to zero for large rt). 


(iii) There exists one purely sinusoidal particular solution, to which any solution 
converges for large t. 


Now the motion of a system of vibrating particles under an external force may be 
modelled by a system of differential equations, of the form 


A,® + A,X + Asx = b, coswt + bz sin wt. (3) 
So analysis for this system parallel to that above for the single equation (2) would 
be very informative. To achieve this we need three things. 


(i) To know that the general solution of (3) is made up of ‘particular solution 
plus complementary function’. 


(ii) To have a condition on the matrices A,, A) and A; which ensures that the 
complementary function is a transient. 


(iii) To know that there is a purely sinusoidal solution of (3), and how to 
calculate it. 


Each of these can be dealt with, and we will do so now. Points (i) and (ii) are 
resolved by Theorems 1 and 2 below. These are stated without proof, although a 
proof of Theorem 1 would be quite straightforward. The proof of Theorem 2 
depends on methods in linear algebra beyond the scope of this course. 


Theorem 1 
Suppose that x, is a particular solution of the system 
Ai + Aok + Agx = h(t). (4) 


Then the general solution of (4) is 
X=X,+X. 


where x, (the complementary function) is the general solution of 
the associated homogeneous system 


Aik + Aak + Agx =0. 


We shall study such systems 
in Unit 24. 
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Theorem 2 
Suppose that each of the matrices A,, A> and A; has constant 
real entries, is symmetric, and has the property that each of its A symmetric matrix A is one 
= : ies G for which AT = A (see Unit 20, 
eigenvalues is positive. Then any solution of the system Su ion 4,4). 
A,X + A.x + Asx =0 


has the property that x tends to 0 as t becomes large. (That is, 
every component of the vector x becomes small as 1 becomes 
large.) 


Exercise 1 


Which of the systems below satisfies the criteria given in Theorem 2 for the solution to be a 
transient? 


@) 43%) + 2%. + 2x, —2 = 
38 + 2%, + 3%, — x1 + 2x, =0 

(ii) -¥, + 3%, + 2%, + 2x, =0 
2%) + 2%, + 3x, =0 

(ili), +%, = 2x) +25 
Ry +X. =x, + 2x. 

[Solution on p, 36] 


Exercise 2 

A particular solution of 
Xp=x.+t 
Xe=x +l 


is not difficult to spot. Do so, and hence find the general solution of the system. (You may 
reuse the solution of Exercise 3(i) of Section 3 here.) 


[Soluron on p. 36} 


4.2 Finding a sinusoidal particular solution 


In (iii) above, we wished to know a purely sinusoidal solution of a system of 
equations of the form (3). In this subsection we shall see how such a solution can 
be calculated. The method uses phasors and is an extension of the procedure in 
Unit 6 for finding a sinusoidal particular solution of a single second-order 
differential equation. We look for a solution of the form x = Re(ze) where z is a 
constant vector with complex entries. Let us see an example. 


Example 2 

Find a sinusoidal particular solution of the system 
¥, +X, =sin2r 
¥2 +X, + 2x, = cos 21. 

Solution 

Step 1: First we write the right-hand sides in phasor form: Phasors were defined in Unit 5 
5, +3) = Re(—ie%¥) and used in Units 6 and 8. 
¥2 + %, + 2x2 = Re(e*"). 5) 

Step 2: Next we look for a solution of the form 
x = Re(ze”"). 

Then ifx=[x, x2)’ andz=[z, =)", we have 
x1 = Re(z,e*"), 

so %, = Re(2iz,e”"), 

and ¥, = Re(—4z,e7"); 


and similarly for x3. 
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Thus to satisfy (5) we need 


—42, + 2iz2 = —i (6) 
42, + 2iz; + 222 =1. (7) 
Equation (7) may be rewritten as 
2izy — 222 = 1. (8) 


Step 3: Thus we have a pair of simultaneous equations, which may be solved by 
Gaussian elimination, as follows. Adding 2i/4 times (6) to (8) gives 


ie. z= +}. 
Substituting in (6) gives 

—42, -i=-i 
ie. 2,=0. 


Step 4: Now z, and z2 are the phasors of the required particular solution. The 
solution is therefore 


x,=0 
x2 = —4c0s 21. 
The method used in the example is set out in the following procedure. 


Procedure 4.2 
To find a sinusoidal particular solution of 
Ai + Ax + A3x = b, cos@t + bz sin wrt. 
1. Write the system in the form 
A,X + Aax + Ax = Re(de). 


2. Look for a solution of the form 
x = Re(ze), 


wherez=([z, 22... 2, and each z, is a complex 
number. This leads to a system of simultaneous linear 
equations for the z,’s. 


3. Solve this system of linear equations, and so find z. 
4. The required solution is then 
x = Re(ze’’). 


The only significant difference between this procedure and that in Unit 6 for a 
single equation is that here, in Step 3, we must solve a set of simultaneous linear 
equations, rather than just a single equation. 


Let us see what happens if we apply Procedure 4.2 in general. Suppose we have 
written our system in the form 

A,X + Ajk + A3x = Re(de) 
where d is some constant, complex, vector. We look for a solution of the form 
x = Re(ze), Then x = Re(iwze™) and ¥ = Re(—w*ze‘). So we have a solution 
so long as 

A, (—w7ze) + Az (ieze') + A3(ze) = de, 
ie. (—w?A, + iwA, + A3jz=d. 
This represents a set of n simultaneous equations (with complex coefficients) in n 
unknowns (compare with Equations (6) and (7) in Example 2). These equations 
will usually have a unique solution for z, and so we will obtain a sinusoidal 
particular solution x = Re(ze'™). 
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There could be problems with this method if the matrix —w?A, + imA, + A; is 
singular. (This compares with the problem of finding a particular solution of, for 
example, the single equation 


¥ + 4x = sin 2r, 


which has no solution of the form x = Re(ze").) So, once again, we have a 
procedure that will usually work, but which could break down. I will not get into 
a discussion of how to deal with possible ‘problem cases’ here. 


The question below involves the use of Procedure 4.2. You may wish to treat it as 
an Exercise, or as another worked Example, before going on to try Exercise 3 for 
yourself. 
Question 
Find a sinusoidal solution of 

¥,+8X.+x,= 8sin31 — 8cos3r 

Ket X = —6sin 31 — 3cos 3t. 
Solution 


Step 1: First we write the right-hand side in the appropriate form 
%, + 8X2 + x, = Re(—B8ie™ — 8e"") = Re ((—8i — 8)e™) 


SH. +%) = Re ((6i — 3)e*"), 
Step 2: Next we look for a solution of the form x = Re(ze*"), where 
x=[x;  x2J’ andz=[z, zz)’. This is a solution if 

(3i)?2) + 8(3i)22 + 21 = —81-8 


(3i)*z2 + 3izy 
That is, if 
—8z, + 24iz. = —81-8 


Biz, - 92, = 61-3. 


6i — 3. 


Step 3: Multiplying the first equation by 3i/8, and adding this to the second 
equation, gives 


— 182, = 31. 

ie. 22 = -ti. 

Substituting this into the first equation gives 
—82,+4= -81-8. 

ie. z=$4i. 


Step 4: Now 2, and z, are the phasors of the required particular solution. The 
particular solution is therefore 

xX, = 300s 3 — sin 3r 

x2 =$sin 31. 


(As usual, this solution should be checked by substitution in the original 
equations.) 


Exercise 3 
Find a sinusoidal particular solution of each of the systems below. 
(i), + 4x, + 2x2 = 6cos 2r 
Xp + x, + 9x2 = 2sin2r 
(i) + 2%) + %:+ x1, —3x,=sint 
3X, + X2 + 2%. +2x, + x2 =cost—2sint 
[Solution on p. 36] 
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Summary of Section 4 
We discussed results about the solution of the system 
A,X + A2x& + Asx = b, coswt + bz sinwt (3) 


relevant to the motion of a system of vibrating particles. The results are similar to 
those for a single second-order equation. The general solution of (3) is the sum of 
any one particular solution and the complementary function: that is, the general 
solution of A, + A2x + A3x =0. We saw that if each of the matrices A;, A, and 
A; has constant real entries, is symmetric, and has all its eigenvectors positive, 
then the complementary function is a transient; that is, all the entries in the vector 
x tend to zero for large t. We saw also that (3) has a particular solution of the 
form x = d, cos wt + d, sin wt. This can be calculated by looking for a solution of 
(3) of the form x = Re(ze') and calculating the complex vector z from the 
resulting linear equations, as described in Procedure 4.2. 


Thus, so long as the matrices A,, A and A; are of the appropriate form, any 
solution of (3) settles down to the same steady-state solution in the long-term. This 
is independent of the initial conditions, which only affect the transient term. 


5 End of unit exercises and problems 


Exercises 1-3 below cover the main objectives of the unit. They can be used either 
for further practice, or when revising as a check that you have covered the central 
ideas of the unit. 


Problems 1 and 2 take the ideas in the unit a little further. 


Exercise 1 
Find the general solution of each of the systems below. 
(i) 2%, = 8x, — 5x2 (vy) -¥, = 8x, — 5x2 + sin2r 
x2 = 10x; — 7x2 ¥; = 10x, — 7x; + 20082 
(ii) 2X, = 8x, — Sxz + 2e (vi) X= —x, + 2x2 
%2 = 10x, — 7x2 — e* %) = —x, — 3x; 
(ili) X, = 6x, — 4x, 3 = x2 — 4x3 
%, = 10x, — 6x, (vii) %, = 4x, + 4x, —%, 
(iv) %, = 8x, — 5x2 Sey + Xg = — 35 tBa 
¥, = 10x, — 7x2 %y— 4xy = 0. 


[Solution on p. 36} 


Exercise 2 
Describe the behaviour of the solutions of the following system for large values of r. 
4X, + ¥, + 3%, + 3x, + 4x, = cos 2r 
Ey + 4X2 + 3X. + 4x, + 6x, = sin 21, 
[Solution on p. 39} 


Exercise 3 


Construct a set of recurrence relations, based on Euler’s method, suitable for generating 
approximations to the solution of the following system. 


X= x2 


X2=X3 

Xs =x, + 3x2 — 0x3 — cost, 
where x,(0)=2, x2(0)=—1, x3(0)=4. 
[Solution on p. 39] 
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Problem 1 

(i) Suppose that y = x, and that x, satisfies the system given in Exercise 3. Show that y 
satisfies the differential equation 

¥+ oH —3y—y= —cost, 
(ii) Can the process in (i) be reversed? Given that y satisfies the differential equation 
¥ + 4ty = sin 2, 
with y(0) = 1, (0) = 2, (0) = 0, can you construct a first-order system that y must 
satisfy? 
Use this idea to construct a set of recurrence relations to generate numerical solutions 
to this differential equation. 

(iii) The idea introduced in (ji) above can be extended to systems of higher-order equations. 
Consider, for example, the pair of second-order (non-linear) equations below, which 
model the motion of a particle in two dimensions under the influence of gravity, and of 
a resisting force whose size is proportional to the square of the velocity. 


aoe = A? + yx 


—ky 
vey 
Here k and 4 are constants. 


By introducing new variables x; = x, x; = y, x3 = 
system equivalent to the second-order system, 


Hence construct a set of recurrence relations that would generate numerical solutions 
to this second-order system (for given values of k and 4). 


[Solution on p. 40] 


J A PY, 


‘, X4 = J, construct a first-order 


Problem 2 (This introduces methods that cover the special cases not dealt with by the 
Procedures in Sections 1 and 2.) 


(i) Show that x = ae“ + bre* is a solution of the system x = Bx so long as 
(a)b is an eigenvector of B corresponding to the eigenvalue 4, and 
(b)a satisfies (B — Al)a = b, 
(ii) (a) Try to solve the system 
X= x2 
Xp = —4x, + 4x 
by Procedure 1.2 of Section 1, What goes wrong? 
(b) Use part (i) to find two linearly independent solutions of this system. 


(c) Write down the general solution of the system. (You may assume that it is a linear 
combination of the two solutions you have found.) 


(iii) Can you find a solution of the form specified in (i) for the system 
x) = Xz 
Xz = —2x, + 3x2? 

(iv) Solve the system 


X= —4x, +4 +e, 


using the steps outlined below. 
(a)In part (ii) you calculated an eigenvector b of the matrix 


and a vector a such that (B — /l)a = b. Form a matrix P whose columns are b and 
a (in that order). Show that 


ml Pina 
pnp =|" ah 


(b)Make the substitution x = Py in the given system, to obtain 
dr = 21 + yr 
Jr =2y2 +e. 
Solve the second of these equations for y, and hence solve the first equation. 
(c) Hence solve the given system. 
[Solution on p. 40] 
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Appendix 1: Solutions to the exercises 


Solutions to the exercises in Section 1 


1. (i) Differentiating the components of 
g(t)= [207 4 — e**}” we obtain 
a) =[4 0 2c] 
(ii) Let a= [a, a, a, )". Then 
x=ae"=[aje“ ae* (rie 
Thus S 
%=[Aaye” —haye* Jae*! 
adel a, ». Gal” 


= jae”. 
(This is just the result we would have hoped for.) 
2. Systems (1), (2) and (3) are linear, (4) is not. None of (1), 
(2) and (3) are homogeneous. (1) and (2) are constant- 


coefficient, (3) is not. (1) and (3) are first-order; (2) is second- 
order, 
3. (i) It is helpful to rewrite the system as 
X, — 3X2 = -2x,-x, 
Xy= Xp $x FH Xy 
=%) 4 %—35— 0: 
In matrix form this is 
Ax=Bx — (or Ax — Bx = 0) 


1-3 0 —2 -1 0 


withA=| 0 0 1 and B=] 1 i 1}. 


feel | » =1 0 0 0 
(ii) It is constant-coefficient and homogeneous. It is first- 
order. 
4. (i) This system is A,X + A2x = h(t) where 


ig es 


We need Aj '. Since A, is 2 x 2, we can just write this down 
{see Unit 20, Subsection 4.1): 


eof 


In normal form the system is 
&= —Aj'A2x + Ay h(t), 


“melt HE Lg 


and 
artmo=[3°* 


So the system, in normal form is 
-1 -4 +e 
ce me 0 

(ii) This system is A,X + A>x = h(r) where 


wef ab sft a} w-[t} 


The rows of A, are linearly dependent, so A, is singular and 
Aj! does not exist. Hence this system cannot be written in 
normal form. 


Ay'A, = 


(iii) This system is not linear, and hence cannot even be 
written in matrix terms, let alone normal form. 


(iv) This is A\X = Aox where 


wf] af 


Now 


so 
e cig ME 
Ay'Ap = 
Pia am =i 

In normal form the system is therefore 


& =A; 'Aox 
where Aj ‘Ao is as above. 
5. The system is already in the normal form x = Bx where 


2 3 

B= i} 
‘We need to find the eigenvalues and eigenvectors of B, 
The eigenvalues are found by solving the characteristic 
equation det(B — 41) = 0. Now 

“ 2—-A 3 

deu(B — at) =| 2 | 
=(2-A)(1-4)-6 
P-3-4 
=(A-4)(A+ 1). 
So the eigenvalues are 4 and —1. 
To find the eigenvectors we solve 


[2 2alll-Le) 


putting 4 equal to each eigenvalue in turn. 


Case 2 = 4 gives 
—lu+3v= 
2u — 3v=0. 
a an ened corresponding to the eigenvalue 4 is 
Case 4 = —1 gives 
3u+3u=0 
2u + 2v=0, 


So an Cigenvector corresponding to the eigenvalue —1 is 
fh =1f. 


We have two linearly independent eigenvectors, as required, 
Hence the general solution is 


ses tp 


xy = 3Cye" + Cre 

X2 = 2Ce* — Cre" 
where C, and C, are arbitrary constants. 
6. (i) We have x = Bx with 


[i a 
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Now 
dea — ay = [> ‘i 
=) =J 
=P -S514+4 
=(4- 1-4). 


So the eigenvalues of B are 1 and 4. 
To find the corresponding eigenvectors we solve 
[a AE -[] 
-1 -ajle}*Lo 
putting 4 equal to each eigenvalue in turn. 


Case 4 = | gives u + v = 0 so an eigenvector corresponding 
to the eigenvalue | is [1 —1]". 


Case 4 = 4 gives u + 4v = 0 so an eigenvector corresponding 
to the eigenvalue 4 is [—4 1]. 


There are sufficient eigenvectors, so the general solution is 
1 —4 
x-o[_jf+ef tk 


x1 = Ce — 4C,e* 
x2 = —Cye + Cye*. 
(ii) We have x = Bx with 


or 


5 -6 -6 
B= |-1 4 2 
3-6 -4 
Now 
S-iA -6 -6 
det(B-Al)=|-1 4-2 2 
$e H¥-3 


= (5— A) [(4 — 4)(—4 — 4) +12) 
+6[4+ 4-6] -6[6 - 12+ 34] 

= -B +52 -B +4 
=-(A-1)4-27. 

So the eigenvalues are 1 and 2 (repeated). 

To find the corresponding eigenvectors we solve 

[5-4 -6 —6 jfu] fo 

-1 4-4 2 v}=|0 

3 -6 -4-A])|w. 0 


putting A equal to each eigenvalue in turn. 


Case A = 1 gives 
[4 -6 -6]fu] fo 
-1 3. 2ilel=lo 
L3 -6 -sllw] Lo 


Performing row operations we obtain 


fo 6 2] [‘ 0 


-1 3 2i|el=lo 
Lo 3 alt Lo 
ro oo fu] fo 
-1 0 1}jv }=|0 
o 3 1flw] Lo 


Therefore 
—u+w=0 
jv+w=0. 


So an eigenvector corresponding to the eigenvalue | is 


[-3 1 -3}7 
Case i = 2 gives 
3-6 -6][u] [fo 
-1 2 2|/v]/=/o} 
3-6 -6}|w} Lo 


This reduces to the single equation 
—ut+2v+2w=0, 


Hence there are two linearly independent eigenvectors, for 
example (2 0 = ‘1? and [2 1 0). 


We have sufficient linearly independent eigenvectors, so the 
general solution is 


=3 2 2 
x=C,| Lle+C,]0}e* + Cy] 1 Je 
-3 i) 0 


or 
X= —3C\e + 2Cze + 2Cye* 
x2 = Cie + Cye™ 
Xs = —3Cye + Cre”. 

(iii) We have x = Bx where 


det(B — AI) = 


2-4 | 
-1 -] 
= +242 
So there is a conjugate pair of eigenvalues —1 + i, 


The equation for the eigenvector corresponding to the 
eigenvalue —1 + i is 


ee 1 * l]-[e} 


(-1—i)u+20=0 
—u+(l—ip=0. 
The first equation is (1 + i) times the second and so both 
give 2v = (1 + i)u. Thus an eigenvector corresponding to the 
eigenvalue A= —1 + iis a= (2 1+ i)". (The eigenvector 
corresponding to the eigenvalue 1 = —1 — i is therefore 
@=(2 1-—i)", but this need not concern us here.) 


Now 


“Lik 
Iti 


mal 2(cost + isint) | 


or 


(1 + i) (cost + isint) 


| 2cost ad 2sint 
cost —sint cost + sint |}’ 


Resa) = ef pcos \ 
cost — sint 
Imac) =f ae } 


cost + sint 


Thus, by Procedure 1.2(a), we can write down the general 
solution in the form 


xe fa 2eost Jee sine ]}. 
cost —sint cost + sinr 
that is © 

x, =e-"(2C, cost + 2C)sint) 

xX, =e-"((C, + Cz)cost + (Cz — C,)sint). 


Solutions to the exercises in Section 2 

1. Note that a good deal of the calculation has already 
been done in Example 1. We have x = Bx + h(t) where 
B= 2 : as in the Example, but now 

hie) = [> 4e**]7. 

Hence Step | of Procedure 2.1 is already done. We have 


J} eel 
4 


We next put x = Py, to obtain 


y=P"'BPy + P~'h(t), 
that is 


Cl-b ICG 12] 


di=4y + e* 

Ya= —Ya— 22%. 
Each of these equations can be solved using the integrating 
factor method. Writing the first equation in the form 

dr —4y, =e, 


we see that the integrating factor is e~“, Multiplying both 
sides by this gives 


oo 


or 


fen) seo “e% met 
so ety, = feta —te* +0, 
ie. Vi = Cie“ — fe, 
Writing the second equation in the form 
Jat Ya = —20%, 
we see that the integrating factor is e’. Multiplying both sides 
by this gives 
sen) = —2e%, 
so 
ey, = —te* + C2, 


y2 = Cze™' — Ge, 
The solution of the original system is x = Py; that is, 


iy -des--tah 


1 = 3Cye" + Cze* — Be” 
X2 = 2Cye" — Cze~t — 4e*. 


or 
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2. We have 
% = Bx + h(t) 
—2 10 
10 -2 
Solving det(B — A) = 0 gives the eigenvalues of B. Now 
i =A 10 
10 -2-4 
=(2+ 47 - 10? 
= (A + 12)(A— 8). 
So the eigenvalues are —12 and 8. 
Eigenvectors corresponding to the eigenvalue 8 satisfy. 
[we -oll:)-[0} 
10 —10jLe 0 


that is u =v. One such eigenvector is [1 1)’. 


where B= | Jana h(t)= [48 91e*]. 


det(B — 


Eigenvectors corresponding to the eigenvalue —12 satisfy 


10 -10][u]_fo 
10 10][v} Lo 
that is w + » = 0. One such eigenvector is [1 


ae | 
1-17 


-1)}". 


tar =[ 


Putting x = Py in the original system gives 
y=P°'BPy + P~'h(t). That is 


Cl-Lo -allI+ lle} 


Di =8y, +244 Ye 
Ja = —12y, + 24 — Be, 


We can solve these equations by the integrating factor 
method. We obtain (omitting the details) 


ys = Cye™ — (3 + Hel) 
ya = Cze“!™ + (2 - Je’). 


_fi tf Cre*-6 + ¥e) 
x=Pr-[} By Pence 


x, = Cye™ + Cze™!* — 1 — 102 

X2 = Cye™ — Ce! — 5 — 3. 
3. (i) As in Example 2, we can write down recurrence 
relations corresponding to Euler's method directly. We 
obtain 

Xie =X 1, + h(3rhX2, + 4) 

Xare1 = Xa, t+ W(hX 1, — X32, — €"). 


These are already in a form suitable for generating solutions, 
although it is a little neater to rewrite the second equation: 


Xaye, =rh?X,, + (1—)X,, — her. 
We have also that Xj.9 = 5, X2,9 =2. 


Gi) With h = 0.1 we need to calculate X;,. and X2,2 to 
approximate x, (0.2) and x2(0.2). Putting r= 0 in the 


or 


or 
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recurrence, we obtain 

Xi, =5+01x4=54 

Xo, =0+09 x 2-01 =17. 
Putting r = 1, we obtain 

X12 = X11. +01(0.3X2, + 4) 
4+0.1 x 451 
= 5.851, 


= 0.01 x 5.4+0.9 x 1.7—O.1e%! 
= 1.473. 


Thus we obtain the approximations x, (0.2) ~ 5.851 and 
X2(0.2) = 1.473. 


Solutions to the exercises in Section 3 
1. To calculate the eigenvalues of A, we first form its 
characteristic equation: 
—$-4 § | 
$  -§-4 
=(4 +87 - 67 
=(A+ 1IA+4), 
giving the eigenvalues —1, —4. 
To find the corresponding eigenvectors we solve 


[3% 3-all)-[] 
$ -$-Allel” Lo 
putting A equal to each eigenvalue in turn, 
Case 4 = —1 gives the equations 
gu + $v = 0 (twice). 
Thus u= —vand so [1 —1)" is an eigenvector 
corresponding to the eigenvalue —1. 
Case 4 = —4 gives the equations 
—gu+gv=0 
gu — $v =0. 
Thus u =v and so [1 1)’ is an eigenvector 
corresponding to the eigenvalue —$. 


2. The pair of differential equations can be written # = Ar 
where 


LY 


Let us first find the eigenvalues and eigenvectors of A. To 
find the eigenvalues, form the characteristic equation. 


0= 


vena od, 
=A+PA+#-# 
= + RA + 
=P +544 


=(4+ 440). 
Thus the eigenvalues of A are —1 and —4. 
The eigenvectors corresponding to —1 are given by 
[37 elle)-Lo} 
$ —¥+1}lel” Lo 
That is 


— Bu + $0 =0 
fu —Fv=0. 


Thus v = 3u and so an eigenvector corresponding to the 
eigenvalue —lis [1 3)”. 


The eigenvectors for —4 are given by 


[3 9 a]-[] 


That is 
jut $v=0 
fu +o =0. 


Thus u = —2v and so an eigenvector corresponding to the 
eigenvalue —4is[-2 1 ]". 

Since the eigenvalues are negative, the object will describe 
simple harmonic motion along the corresonding eigenvectors. 
The frequencies will be the square roots of minus the 
eigenvalues. So there are two possibilities: 


(i) motion in the direction of the vector [1 3)", with 
angular frequency 1; 
(ii) motion in the direction of the vector [—2 1)", with 


angular frequency 2. 
3. (i) We have ¥ = Bx where 
0 1 
e-[I 
From the characteristic equation 


ee a) a 


The eigenvalues of B are | and —1. 
The eigenvectors corresponding to the eigenvalue 1 satisfy 


= 1|fu 0 
1 -1]Le}~ Lo 
So one such eigenvector is [1 Tins 
The eigenvectors corresponding to the eigenvalue —1 satisfy 


1 1]fu 0 
[i Cl-Leh 
So one such eigenvector is (1 —1]". 
By Procedure 3.2 the general solution is 


x= [ove + Dye") + [ 4 fesoose + Dzsint) 


A 
|-#-1 


or 
xX, = Cyeé + Dye“ + C, cost + Dz sint 
X2 = C,é + Dye“! — C, cost — D, sint. 


(ii) We must first divide through the second equation by 2, 
to obtain = Bx with 


=10 -6 
a-[ -6 ex 
From the characteristic equation 
-10-4 -6 | 
=6 -10-4| 
=(+ 107-6 
= (4 + 16)(A + 4) 
the eigenvalues of B are —4 and —16. 
The eigenvectors corresponding to the eigenvalue —4 satisfy 


[5 =<l[-)-[0] 


So one such eigenvector is [1 —1]". 


o-| 


The eigenvectors corresponding to the eigenvalue — 16 satisfy 
[-< “<ll-}-Lo} 
-6 6] 0 
So one such eigenvector is [1 1]". 
By Procedure 3.2 the general solution is 


x =[_{]evcos2 + D, sin 2t) 


+ [| |cssosae + D,sin4t), 


or 
xX, = C, cos 2t + D, sin 2t + C, cos 4t + Dz sin 4r 
xX, = —C, cos 2t — D, sin 2t + C, cos 4t + D3 sin 4t. 


Solutions to the exercises in Section 4 


1. In each case we must write the equation in the form 
A,X + A2% + Asx = 0, and see whether the matrices A,, A> 
and A, Satisfy the specified criteria. 


(i) Here 
i 0 Se: 
ae [0 } aan [> A} 


Each of these matrices is symmetric. A, has eigenvalues | 
and 3. We must find the eigenvalues of A; and A;. 


For A; we have 


B—A 2 
devia, ~ 2) =| ee 
=(3-Ap-2? 
=(A-1)(A—-5). 


So A, has eigenvalues 1 and 5. 
For A; we have 


@-a 1 
detias — ay) =P ay 
= (2-4? —-(-1) 
=(4-1)—3). 


So As has eigenvalues 1 and 3. 

Hence all the eigenvalues of each matrix are positive. 

So the conditions of Theorem 2 are satisfied in this case. 
(ii) Here 


1 0 er 3 
el CO ce FO 
2 #0 
ali ol 
Each of these matrices is symmetric, and A, and A; have 


positive eigenvalues. However A; has 0 as an eigenvalue, so 
the conditions of Theorem 2 are not satisfied in this case. 


{In fact, one can see that x; = 0, x2 = a (where a is a non- 
zero constant) is a solution of this system—and is not a 
transient. } 

(iii) To compare with the form in Theorem 2 we need to 
bring all the terms to the left-hand side of the equation . So 


here 
—2 -1 
a 
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Each of these matrices is symmetric, and A, and A; have 
positive eigenvalues. However, 


det(As — Al) = (—2 — 4)? — (-1)7 
=(-1-4)(-3 - 4). 


Hence the eigenvalues of A, are negative (—1 and —3), and 
the conditions of Theorem 2 are not satisfied in this case. 


2. Ifx, = —1, x; = —1, then ¥, = ¥, = 0, and the equations 
are satisfied. 


From Solution 3(i) in Section 3, the general solution of 
#, =x 
y= x 
—the associated homogeneous system—is 
x, = Cye + Dye! + C, cost + Dy sint 
x2 = Ce + Dye“ — C, cost — Dy sint. 
So the general solution of the given system is 
x, = Cye + Dye + C, cost + Dz sint,— 1 
X2 = Cie + Dye™' — C, cost — Dz sint — t. 
3. (i) We have a solution of the form x = Re(ze*") if 
—42z, + 42, +22, =6 
— 422 + 2) + 922 = —2i. 
That is, if 
222 =6 
2 + Szz = —2i. 


Thus z; = 3, z, = —2i — 15 and so the required particular 
solution is 

x, = 2sin 2t — 15cos 21, 

X2 = 3cos 21. 


(ii) We have a solution x = Re(ze*) if 
—2, — 22, + iz, + 2) — 32. = -i 
—32, — 22 + 2izz + 22, + 22 = 14+ 21 

That is, if 
iz, —5z,3=-i 
—z, + 2iz, = 1+ 21. 


Adding i times the second equation to the first gives 
—7z; = —2,s0 


n=} 
From the second equation 
n= ¢ (142i 
=-1--—- 
The required particular solution is therefore 


x, = —cost + 4fsinr, 
X2 = Foose. 


Solutions to the exercises in Section 5 


1.3) The system can be written x = Bx with 


e-[ 5 7 


To find the eigenvalues of B we solve the characteristic 
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equation det(B — AI) = 0. We have 
Is—A =5 | 
10. =T=2 
2-2-6 
=(A-—3)(A+ 2), 
80 the eigenvalues of B are 3 and —2. 
The eigenvectors corresponding to the eigenvalue 3 satisfy 
[iw =ll:]-[c} 
10 —10}| 0 0 
Thus u = v and so an eigenvector is [1 1h ea 
The eigenvectors corresponding to the eigenvalue —2 satisfy 
[wm =sI[-)-[o] 
10 -—SjLv 0 
Thus » = 2u and so an eigenvector is [1 2)”. 
By Procedure 1.2 the general solution of x = Bx is 


x= oft feels} 


That is 
x, = Cye™ + Cye"** 
x2 = Cye™ + 2C,e~. 


dexn — a) =| 


(ii) This can be written x = Bx + h(t) where B is the same 
as in (i) and h(t) = (2e' —e**]’. The general solution can 
be found using Procedure 2.1. 


The eigenvalues and eigenvectors of B are given in (i). 
We form the matrix 
1 1 
P= 
[4] 


and put x = Py. This gives 


y= P-inry spf 

fs OL.f 2 -1f 2 

“10 =9F Thad alh—af 
3y1 + 4e + 2 
—2y, — 2e — e*, 


Using the integrating factor method we find that the-general 
solutions of these equations are 


yy = Cye™ — 2e — e* 

y2 = Cye"™ — Fe — fe”. 
So the general solution of the original system is x = Py, that 
Is 

x, = Cye™ + Cye-™ — $e — je* 

Xz = Cye™ + 2C,e~4 — Wet — ge". 


(iii) This can be written & = Bx where 


e-[15 Wel 


The eigenvalues and eigenvectors of this matrix are complex, 
One eigenvalue is A = 2i, with corresponding eigenvector 
a=[4 6-2i]". 
Now 
os 4 |e 4(cos 21 + isin 2t) 
ae! = ler = z 
6—2i (6 — 2i) (cos 2r + isin 2r) 


so 
4008 2t 
Resse) =| a ee 
4sin2t 
tina = Vée 2r- Sadek 


By Procedure 1.2(a) the general solution is 

x = C,Re(ae“) + CzIm(ae*), 
or 

%, = 4C, cos 2t + 4C, sin 21 

Xz = C,(6cos 2t + 2sin 2t) + C,(6sin 2t — 2cos 2t), 
(iv) In matrix form this can be written % = Bx where B is 


the same matrix as in (i). The eigenvalues and eigenvectors of 
B are given in (i), By Procedure 3.2 the general solution is 


x= [i Joe + Dye“) 


+ [5 ]cssn /2t + Dz cos ,/21) 
or 
x1 = Cre" + Dye" + Cz sin /2t + D3 cos /2t. 
x2 = Cie" + Dye“ + 2C sin \/2t + 2D; c08 \/21. 


(vy) The general solution is the sum of the complementary 
function, which we found in solution (iv), and any one 
Particular solution. We can find a sinusoidal particular 
solution by applying Procedure 4.2. 


Putting x = Re(ze*") leads to the equations 
—42, = 82, — $2, -i 
—4z; = 102, — 723 +2, 
which have the solution 
n=-$-h, = —H- 91 
Hence a particular solution of the system is 
—4oos 2t + yysin 2 
—Yoos2r + $sin 21. 
The general solution of the system is therefore 
x = Cyev + Dye™*™ + Cy sin /2t + Dz cos /2t 
— $c0s 2r + fxsin 2¢ 
x2 = Cie + Dye~¥™ + 2C; sin /2t + 2D, 008 ,/2t 
— Y¥cos2r + 4sin 2r. 
(vi) This equation can be written x = Bx where 


x 


*2 


=k) 2) 8) 
B=/-1 -3 0}. 
0 1 -4 


To find the eigenvalues of B we solve the characteristic 
equation det(B — Al) = 0. Now 


-A 2 0 
de(B-A=| -1 -3-4 0 


=(-4-4) 


=i -3-j 
=(-4-A(#? +4045). 
So the eigenvalues of B are —4 and —2 +i. 
The eigenvectors corresponding to the eigenvalue —4 satisfy 
3 2 Off 0 
au 1 O}| v |=] 0}. 
0 1 O}Lw. 0 


So u=v =O, and an eigenvector is (0 0 1). 


The eigenvectors corresponding to the eigenvalue —2 + i 
satisfy 


1-i ul] [o 
-1 v}| =| 0 
0 w} LO 
That is 
(1—iu+2v=0 
—u-(1+ijp=0 
v—(2+iw 


The first equation is (i — 1) times the second equation. The 
second equation gives u = —(1 + i)v and the third equation 
gives w = v/(2 + i) = 4(2— ie. Thus an eigenvector 
corresponding to the eigenvalue 4 = —2 + iis 
a=(-(i+) 1 4(2-H]" 
The general solution of the given system can be found using 
Procedure 1.2(a). We have 
ba 
acme 1 
42-3) 


(cos t + isint) 


—cost + sint 
cost +i sint 


—sinf — cost 
se* 
4(2cost + sint) 


4(2sinr — cost) 


so 


cost + sin 
Re(ae") =e" cost q 


4(2cost + sin 


—sint — cost 
| sint 


4(2sint — cost) 


Im(ae“) = e~ 


The general solution is therefore 


0 —cost +sint 
x=C,/0le“+C cost ent 
1 4(2cost + sint) 
—sint — cost 
+C3 sint ex, 


4(2sint — cost) 


or 
x1 = C,e72"(—cost + sin t) + Cye~?"(—sint — cost) 
Xz = C2e7* cost + Cse™*sint 
x3 = Cie + 4C3e°*(2cost + sint) 
+ 4Cye~*(2sint — cost). 
(vii) We first write the equations in standard form. To do 
this, we separate the derivative and non-derivative terms: 
My +%) = 4x, + 4x2 
y= 3x, + Xp+ Xs 
X= 4x3. 
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and then write this in the matrix form Ax = Cx where 


t 6 4 4 «6 
A=|0 1 O| and C=/3 1 1 
oo 1 0 0 4 


The norma! form of the equation is 
*%=A Ce 


The inverse A~' can be found using row operations (see Unit 
20, Subsection 4.2). We have 


A 1 
1 1 0 1 0 O7R, 
0 1 0 0 1 O}R; 
0 0 1 or 0 Rs 
Ri=ki (6 ~ ola =e 
0 1 0 0 1 0 
0 1 i 1 
I aS 
Thus 
1 a | 0 
A‘=|0 1 0 
o oO 1 


So the normal form is 
1-1 oyf4 4 £0 
x=/0 1 0}|3 1 1\x 


[Oe ee | Ce 
1 3 <2 

=|3 1 1)x. 
0 0 4 


We next need the eigenvalues and eigenvectors of 


1 2 
B=|3 1 1}. 
0 Oo 4 


These were calculated in Example 1 of Section 3 where we 
found 


eigenvalues | corresponding eigenvectors 


By Procedure 1.2 the general solution is 
1 1 0 
x=Cye™ 1-1 14+ Cye“]1] + Cre] 1), 
0 i} 3 


or 
xX; = Cie" + Cre 
x, = —Cye"* + Cre + Cae 
x3 =3C3e". 
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2. The system can be written in the form 
A,X + Ax + Agx =hi(t) 
with 


Each of the matrices A;, Az and A; is symmetric. The 
eigenvalues of A, are 3 and 5 for 

(det(A, — Al) = (4 — AP? — 1? = (A—3) (4-5), 
and those of A; are 3 and 3— all positive, Those of A, are a 
little harder to find. We have 


det(Ay — Al) = (3 — 4)(6 — 4) — 16 
=# -94+2. 


So the eigenvalues of A; are $(9 + ,/9? — 8). Written in this 
form we can see that these eigenvalues are positive, since 
VP -8 <9. 


Hence we may deduce from Theorem 2 of Section 4 that the 
complementary function of the given system is a transient. 


There is a unique sinusoidal particular solution of the given 
system. Every solution of the given system converges to this 
same solution as ¢ becomes large, Let us now calculate this 
solution. 


39 


To do this we put x = Re(ze"). This yields the equations 
—16z, — 42; + 6iz, + 32, +422 =1 
—4z, — 16z, + 6iz2 + 42, + 62, = 
That is 
(61 — 13)z, 


61+ 13 

205 * 
_ _ (61+ 10) —6 + 101 
ae” Tian ia 


The steady-state solution is therefore 

X; = —xbs(13 cos 2r — 6 sin 21) 

X2 = —rhe(6c0s 21 + 10sin 21), 
3. The relevant recurrence relations, for step-length h, may 
be written down directly: 

Xirer =Xi, + hXa, 

Kaper = Xr, + hXy, 

Xayor = Xa, + WX, + 3X2, —rhXy,, — cosh) 
where X;9=2, X29=—-1. X30=4 


(Here X,=[X,, X2,  X3,]' approximates the true 
value of x(rh).) 
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Appendix 2: Solutions to the problems 


Solutions to the problems in Section 5 
1. (i) We have y = x), ¥; = x2, X; = x3, and 
Xy = xX, + 3x2 — tx3 — cost. i) 


We can express the X3, x), x2 and x; in (1) in terms of y as 
follows 


“u=y) 
*2=%1 =), 
x3 =X. = 5, 
=i 


Substituting these into Equation (1) gives 
V=y +39 — 1 — cost, 

or 
V+) — 3) —y = —cost, 

as required. 


(ii) We are given that y satisfies’ + 4ty = sin 2r. To reverse 
the process in (i), we introduce some new variables. We 
define x,, x) and x; by the equations 


ay 
m=) 
X= ip. 
Then we have 
X= 
and 
Ry = xy. 
Also %; = ¥; so from the given differential equation, 
Xs; = sin 2t — 40x,. 
Thus the given equation is equivalent to the first-order 
system: 
X= x2 
Xa = xy 
Xy = sin 2t — 4tx,. 
The initial conditions give 
*(0)=1, x20) =2, x30) = 0. 
We can use the method described in Subsection 2.2 to write 
down recurrence relations to generate approximate solutions 
to this system: 
Xieer =Xig + hX2, 
Kayes = Xap + hXg, 
Xsye1 =X, + A(sin 2rh — 4rhX ,,,). 
with = X49 =1, X29 =2, X30 =0. 
If we use these recurrence relations, then X,,, will provide an 


approximation to the solution y(rh) of the original 
differential equation. 


(The method described in this problem generalises the 
method of numerical approximation for a single second-order 
differential equation described in Unit 6.) 


(ii) If we put x, =x, x,=y, x3=%, xy we have 
X= X3 
Xp =X 
—kx, 
%3(=%) = — (xd + x3)!2x. 
? xi +33 a 
sy - xd + 99). 
xi +x3 


The following recurrence relations are obtained if we apply 
Euler's method to each equation in turn (with step length h). 


Xiper =Xiy + hXs, 
Mayer =X2,+hXa, 


kX. 2 
Xsy01 = Xa [tee 44X53, + X4,1°%,| 
kX 2, ss 
Kaper = Kar gt + 1003, + XEN] 
2 @ If 
x = ae” + bre” (10) 
then 


& = Jae” + ble" + die), 
So (1) is a solution of x = Bx so long as 
(Ja + bye + bite” = Bae” + bre“), 
This is the case so long as 
a + b = Ba (equating coefficients of e”), (2) 
and i4b=Bb (equating coefficients of re”). (3) 


That is, (a) b is an eigenvector of B, with eigenvalue 2 (from 
(3)), and (b) a satisfies (B — A1)a = b (from (2)), as required. 


(ii) (a) We have x = Bx with 


1 
4-] 


Thus 2 is the only eigenvalue, The corresponding 
eigenvectors satisfy 


[= ICI-(] 

-4 2jLe oy 

Thus p = 2u, and so every eigenvector is of the form 

(kK 2kJT=ktL 2)", 

To be able to apply Procedure 1.2 of Section 1 we require 


two linearly independent eigenvectors, but this is impossible 
since all eigenvectors are multiples of [1 Tits 


=P -4444= (4-2). 


(b) Since [1 2)” is an eigenvector corresponding to the 
eigenvalue 2 we know that 


“LF 


is one solution of x = Bx. We can use part (i) of the problem 
to find a second linearly independent solution. To do this we 
need to find a vector a satisfying (B — Al)a = b, where b is 
the eigenvector [1 
need 


2]'. That is, ifa=([a, a3)" we 


These two equations are equivalent, so there is a wide choice 
of solutions. One is a, =0, a,=1 givinga=(0  1)?. 
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With these values for a and b, we know from part (i) that 


[hse 


is a solution of % = Bx. 


(©) The general solution is 
1 0 i We 
-ofjJe+e([le+ [> fe") 


Xy = Cye" + Care” 
Xz = 2Cye" + Cre" + 2Cyte™. 


(To get this we assume the fact that the general solution is a 
linear combination of the two solutions we have found.) 


or 


(ii) We now have x = Ax, with 


0 iT 
oS [ = “| 
The eigenvalues of A satisfy 
A 1 
-2 3-4 
Thus the eigenvalues are 2 and 1. 


O= 


3h +2=(4-2A- 1). 


‘The eigenvectors corresponding to the eigenvalue | satisfy 


[3 1 A 0 
-2 2flo}~Lof 
Thus wu = v, and so an eigenvector is [1 1)". 


Let us try to find a solution of the form specified in part (i). 
We would need a vectora=[a, a)" satisfying 


[= aM)-0) 
-2  2jlas}~ [1] 
That is 
a,—a,=1 
2a; — 2a, = 1. 
These equations are inconsistent and have no solutions. 
If we proceed similarly for the eigenvalue 2 we again arrive 
at a set of inconsistent equations for a. So in this case no 
solution of the form specified in part (i) can be found. This is 


just as well, as we have enough independent solutions 
already! 


(iv) (a) Inpart(ii)wehadb=(1 2].a=(0 177. 
[You may have different vectors here—there was more than 
one valid choice of these vectors.] Thus 


San 
wmel iles all l 


(b) The given system can be written 


x-o+(°] 


Putting x = Py gives 


nm 


so 

y= P'apy+e[ 

é 

ie. 

et 0 

c [0 ib i [-} 
Thus 

Fi =i +2 

yr= ate. 


We can solve the second equation, by the integrating factor 
method to obtain 


yn = Ce" —e. 
Substituting this into the first equation gives 
dy =2y, + Cye* - 


We can now solve this by the integrating factor method, to 
obtain 


Ys = Cre™ + Cyte™ + 
(c) Now x = Py, so we have 


1 Of Cre + Cyte* + 
2 1 Ce" 


Xy = Cre" + Cyte* +e 
X2 = 2Cze™ + Cye*(1 + 2t) +e 
is the general solution of the given system. 


ig 
[ooee eRe 


